﻿#pragma once

#include <cassert>
#include <cfloat>
#include <climits>
#include <cmath>
#include <initializer_list>
#include <iostream>
#include <vector>

namespace CAGD
{

/**
 * @brief 数学向量类型
 *
 */
class Matrix;
class Vector
{
  protected:
    const int m_dim;
    double *m_value;

  public:
    /**
     * @brief 初始化构造函数.
     * @param D
     * @param t
     * @see
     */
    explicit Vector(int D, double t = 0);

    /**
     * @brief 可以通过列表初始化.
     * @param l
     * @see
     */
    Vector(std::initializer_list<double> l);

    /**
     * @brief 拷贝构造.
     * @param rhs
     * @see
     */
    Vector(const Vector &rhs);

    /**
     * @brief 获得向量长度.
     * @return
     * @see
     */
    int Size() const;

    /**
     * @brief 返回引用.
     * @return
     * @param i
     * @see
     */
    double &operator[](int i);

    /**
     * @brief 返回值.
     * @return
     * @param i
     * @see
     */
    double operator[](int i) const;

    /**
     * @brief 加法.
     * @return
     * @param rhs
     * @see
     */
    Vector operator+(const Vector &rhs) const;

    /**
     * @brief 减法.
     * @return
     * @param rhs
     * @see
     */
    Vector operator-(const Vector &rhs) const;

    /**
     * @brief 内积.
     * @return
     * @param rhs
     * @see
     */
    double operator*(const Vector &rhs) const;

    /**
     * @brief 伸缩.
     * @return
     * @param t
     * @see
     */
    Vector operator*(double t) const;

    /**
     * @brief 伸缩.
     * @return
     * @param t
     * @see
     */
    Vector operator/(double t) const;

    /**
     * @brief 乘矩阵.
     * @return
     * @param rhs
     * @see
     */
    Matrix operator*(const Matrix &rhs) const;

    /**
     * @brief 赋值.
     * @param rhs
     * @see
     */
    void operator=(const Vector &rhs);

    /**
     * @brief 比较.
     * @return
     * @param rhs
     * @see
     */
    bool operator==(const Vector &rhs) const;

    /**
     * @brief 取反操作.
     * @return
     * @see
     */
    Vector operator-() const;

    /**
     * @brief 冒泡排序.
     * @see
     */
    void Sort();

    /**
     * @brief 元素求和.
     * @return
     * @see
     */
    double Sum() const;

    /**
     * @brief 取绝对值.
     * @return
     * @see
     */
    Vector Abs() const;

    /**
     * @brief 获取符号向量.
     * @return
     * @see
     */
    Vector Sgn() const;

    /**
     * @brief 向量范数.
     * @return
     * @param nt
     * @see
     */
    double Norm(int nt = 2) const;

    /**
     * @brief 取最小元.
     * @return
     * @see
     */
    double Min() const;

    /**
     * @brief 取最大元.
     * @return
     * @see
     */
    double Max() const;

    /**
     * @brief 取具有最大模的分量.
     * @return
     * @see
     */
    int MaxId() const;

    /**
     * @brief 转化成矩阵.
     * @return
     * @see
     */
    Matrix ToBMatrix() const;

    /**
     * @brief 转置成矩阵.
     * @return
     * @see
     */
    Matrix Transpose() const;

    /**
     * @brief Householder 变换.
     * @return
     * @param x
     * @param v
     * @see
     */
    static double HouseHolder(Vector &x, Vector &v);

    /**
     * @brief Givens 变换.
     * @return
     * @param v
     * @param i
     * @param j
     * @see
     */
    static Vector Givens(Vector &v, int i = 0, int j = 1);

    /**
     * @brief 隐式类型转换.
     * @return
     * @see
     */
    operator std::vector<double>() const
    {
        std::vector<double> r(m_dim);
        for (int d = 0; d < m_dim; d++)
            r[d] = m_value[d];
        return r;
    }

    /**
     * @brief 输出函数.
     * @return
     * @param os
     * @param v
     * @see
     */
    friend std::ostream &operator<<(std::ostream &os, const Vector &v)
    {
        os << "{";
        for (int i = 0; i < v.Size(); i++)
            os << v[i] << " ";
        os << "\b}" << std::endl;
        return os;
    }

    /**
     * @brief 析构.
     * @see
     */
    ~Vector();
};

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * @brief 数学矩阵类型
 *
 */
class Matrix
{
  protected:
    const int m_row, m_col;
    double **m_value;

  public:
    /**
     * @brief 友元.
     *
     */
    friend class LinearSolver;

    /**
     * @brief 初始化构造函数.
     * @param Row
     * @param Col
     * @param t
     * @see
     */
    Matrix(int Row, int Col, double t = 0);

    /**
     * @brief 通过列表初始化.
     * @param Row
     * @param Col
     * @param l
     * @see
     */
    Matrix(int Row, int Col, std::initializer_list<double> l);

    /**
     * @brief 拷贝构造.
     * @param rhs
     * @see
     */
    Matrix(const Matrix &rhs);

    /**
     * @brief 单位阵.
     * @return
     * @param n
     * @see
     */
    static Matrix Identity(int n);

    /**
     * @brief 行交换矩阵.
     * @return
     * @param n
     * @param i
     * @param j
     * @see
     */
    static Matrix ExchangeRC(int n, int i, int j);

    /**
     * @brief 行乘矩阵.
     * @return
     * @param n
     * @param i
     * @param t
     * @see
     */
    static Matrix TimesRC(int n, int i, double t);

    /**
     * @brief 行加乘矩阵.
     * @return
     * @param n
     * @param i
     * @param j
     * @param t
     * @see
     */
    static Matrix TimesAddRC(int n, int i, int j, double t);

    /**
     * @brief 获得行数.
     * @return
     * @see
     */
    int Row() const;

    /**
     * @brief 获得列数.
     * @return
     * @see
     */
    int Col() const;

    /**
     * @brief 返回 i,j 元素.
     * @return
     * @param i
     * @param j
     * @see
     */
    double operator()(int i, int j) const;

    /**
     * @brief 返回 i,j 元素引用.
     * @return
     * @param i
     * @param j
     * @see
     */
    double &operator()(int i, int j);

    /**
     * @brief 加法.
     * @return
     * @param rhs
     * @see
     */
    Matrix operator+(const Matrix &rhs) const;

    /**
     * @brief 减法.
     * @return
     * @param rhs
     * @see
     */
    Matrix operator-(const Matrix &rhs) const;

    /**
     * @brief 乘法.
     * @return
     * @param t
     * @see
     */
    Matrix operator*(double t) const;

    /**
     * @brief 除法.
     * @return
     * @param t
     * @see
     */
    Matrix operator/(double t) const;

    /**
     * @brief 幂运算.
     * @return
     * @param N
     * @see
     */
    Matrix operator^(int N) const;

    /**
     * @brief 乘向量.
     * @return
     * @param v
     * @see
     */
    Vector operator*(const Vector &v) const;

    /**
     * @brief 乘矩阵.
     * @return
     * @param rhs
     * @see
     */
    Matrix operator*(const Matrix &rhs) const;

    /**
     * @brief 直接赋值.
     * @param rhs
     * @see
     */
    void operator=(const Matrix &rhs);

    /**
     * @brief 比较.
     * @return
     * @param rhs
     * @see
     */
    bool operator==(const Matrix &rhs) const;

    /**
     * @brief 取反操作.
     * @return
     * @see
     */
    Matrix operator-() const;

    /**
     * @brief 矩阵范数.
     * @return
     * @param nt
     * @see
     */
    double Norm(int nt = 0) const;

    /**
     * @brief e 范数 / f 范数.
     * @return
     * @param c
     * @see
     */
    double Norm(char c) const;

    /**
     * @brief 计算条件数，基于无穷范数.
     * @return
     * @see
     */
    double Condition() const;

    /**
     * @brief 计算矩阵行列式.
     * @return
     * @see
     */
    double Det() const;

    /**
     * @brief 计算矩阵迹.
     * @return
     * @see
     */
    double Trace() const;

    /**
     * @brief 获得对角向量.
     * @return
     * @see
     */
    Vector Diag() const;

    /**
     * @brief 获得上三角.
     * @return
     * @see
     */
    Matrix Upper() const;

    /**
     * @brief 获得下三角.
     * @return
     * @see
     */
    Matrix Lower() const;

    /**
     * @brief 转置.
     * @return
     * @see
     */
    Matrix Transpose() const;

    /**
     * @brief 矩阵求逆，flag 返回是否可逆.
     * @return
     * @param flag
     * @see
     */
    Matrix Invert(int *flag = nullptr) const;

    /**
     * @brief 列主元 LU 分解，传入下标数组，返回置换后的数组，flag 返回矩阵的秩.
     * @return
     * @param perm
     * @param flag
     * @see
     */
    Matrix LU(int *perm = nullptr, int *flag = nullptr) const;

    /**
     * @brief Cholesky LL'.
     * @return
     * @param flag
     * @see
     */
    Matrix Cholesky(int *flag = nullptr) const;

    /**
     * @brief 优化的 Cholesky LDL'.
     * @return
     * @param flag
     * @see
     */
    Matrix CholeskyD(int *flag = nullptr) const;

    /**
     * @brief QR 分解.
     * @return
     * @param d
     * @see
     */
    Matrix QR(Vector &d) const;

    /**
     * @brief 上 Hessenberg 化，在 Q 中累计正交变换.
     * @return
     * @param Q
     * @see
     */
    Matrix Hessenberg(Matrix *Q = nullptr) const;

    /**
     * @brief 对称阵化为三对角对称阵，在 Q 中累计正交变换.
     * @return
     * @param Q
     * @see
     */
    Matrix Tridiagonal(Matrix *Q = nullptr) const;

    /**
     * @brief 矩阵上二对角化，在 U,V 中累计正交变换.
     * @return
     * @param U
     * @param V
     * @see
     */
    Matrix Didiagonal(Matrix *U = nullptr, Matrix *V = nullptr) const;

    /**
     * @brief 输出函数.
     * @return
     * @param os
     * @param m
     * @see
     */
    friend std::ostream &operator<<(std::ostream &os, const Matrix &m)
    {
        os << "{" << std::endl;
        for (int i = 0; i < m.Row(); i++)
        {
            os << "{";
            for (int j = 0; j < m.Col(); j++)
                os << "\t" << m(i, j);
            os << "\t}" << std::endl;
        }
        os << "}" << std::endl;
        return os;
    }

    /**
     * @brief 释放资源.
     * @see
     */
    ~Matrix();

  protected:
    /**
     * @brief 计算指定范围的 Householder 变换，作用结果为 HA.
     * @param startOfRow
     * @param endOfRow
     * @param v
     * @see
     */
    void HA(int startOfRow, int endOfRow, Vector &v);

    /**
     * @brief 计算指定范围的 Householder 变换，作用结果为 AH.
     * @param startOfCol
     * @param endOfCol
     * @param v
     * @see
     */
    void AH(int startOfCol, int endOfCol, Vector &v);

    /**
     * @brief 计算对称阵指定范围的 Householder 变换，作用结果为 HAH.
     * @param startOfRow
     * @param endOfRow
     * @param v
     * @see
     */
    void HAHs(int startOfRow, int endOfRow, Vector &v);

    /**
     * @brief 计算指定分量的 Givens 变换，得到 GA.
     * @param i
     * @param j
     * @param cs
     * @see
     */
    void GA(int i, int j, Vector &cs);

    /**
     * @brief 计算指定分量的 Givens 变换，得到 AG.
     * @param i
     * @param j
     * @param cs
     * @see
     */
    void AG(int i, int j, Vector &cs);

    /**
     * @brief 计算指定分量的 Givens 变换，得到 AGT.
     * @param i
     * @param j
     * @param cs
     * @see
     */
    void AGT(int i, int j, Vector &cs);

    /**
     * @brief 双重步 QR 迭代，l,m 是 H22 前面和后面分块的阶数.
     * @param l
     * @param m
     * @param Q
     * @see
     */
    void DoubleQR(int l = 0, int m = 0, Matrix *Q = nullptr);

    /**
     * @brief Wilkinson QR 迭代，针对三对角对称矩阵.
     * @param l
     * @param m
     * @param Q
     * @see
     */
    void WilkinsonQR(int l = 0, int m = 0, Matrix *Q = nullptr);

    /**
     * @brief Wilkinson 位移的 SVD 迭代.
     * @param l
     * @param m
     * @param U
     * @param V
     * @see
     */
    void WilkinsonSVD(int l = 0, int m = 0, Matrix *U = nullptr, Matrix *V = nullptr);

    /**
     * @brief 对称矩阵绕 p,q 分量旋转 JT*A*J，Q 累计 J.
     * @param Q
     * @param p
     * @param q
     * @see
     */
    void Rotate(Matrix &Q, int p, int q);
};

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * @brief 稀疏矩阵类型
 *
 */
class Sparse
{
  protected:
    const int m_col; /// 固定列数
    const int m_row; /// 固定行数
    int *m_IA;       /// 每个元素记录每行开始位置的全局索引，最后一个元素记录有效元个数

    std::vector<int> m_JA;    /// 每个元素记录列标
    std::vector<double> m_AA; /// 每个元素记录矩阵元素

  public:
    /**
     * @brief 构造.
     * @param row
     * @param col
     * @see
     */
    Sparse(int row, int col);

    /**
     * @brief 拷贝构造.
     * @param spm
     * @see
     */
    Sparse(const Sparse &spm);

    /**
     * @brief 获得矩阵列数.
     * @return
     * @see
     */
    int Col() const;

    /**
     * @brief 获得矩阵的行数.
     * @return
     * @see
     */
    int Row() const;

    /**
     * @brief 获得有效元素的个数.
     * @return
     * @see
     */
    int Nnz() const;

    /**
     * @brief 获得存在的元素，不存在的元素认为是零.
     * @return
     * @param i
     * @param j
     * @see
     */
    double operator()(int i, int j) const;

    /**
     * @brief 在指定位置插入元素.
     * @param i
     * @param j
     * @param value
     * @see
     */
    void Set(int i, int j, double value);

    /**
     * @brief 在指定位置加.
     * @param i
     * @param j
     * @param value
     * @see
     */
    void Add(int i, int j, double value);

    /**
     * @brief 拷贝赋值.
     * @param spm
     * @see
     */
    void operator=(const Sparse &spm);

    /**
     * @brief 加法.
     * @return
     * @param spm
     * @see
     */
    Sparse operator+(const Sparse &spm) const;

    /**
     * @brief 减法.
     * @return
     * @param spm
     * @see
     */
    Sparse operator-(const Sparse &spm) const;

    /**
     * @brief 乘法.
     * @return
     * @param k
     * @see
     */
    Sparse operator*(double k) const;

    /**
     * @brief 除法.
     * @return
     * @param k
     * @see
     */
    Sparse operator/(double k) const;

    /**
     * @brief 稀疏矩阵乘向量.
     * @return
     * @param v
     * @see
     */
    Vector operator*(const Vector &v) const;

    /**
     * @brief 稀疏矩阵乘稀疏矩阵.
     * @return
     * @param spm
     * @see
     */
    Sparse operator*(const Sparse &spm) const;

    /**
     * @brief 转换为一般矩阵.
     * @return
     * @see
     */
    Matrix ToBMatrix() const;

    /**
     * @brief 转置矩阵.
     * @return
     * @see
     */
    Sparse Transpose() const;

    /**
     * @brief 释放资源.
     * @see
     */
    ~Sparse();

    /**
     * @brief 输出函数.
     * @return
     * @param os
     * @param m
     * @see
     */
    friend std::ostream &operator<<(std::ostream &os, const Sparse &m)
    {
        os << "{" << std::endl;
        for (int i = 0; i < m.Row(); i++)
            for (int k = m.m_IA[i]; k < m.m_IA[i + 1]; k++)
                os << "(" << i << ", " << m.m_JA[k] << ") -> " << m.m_AA[k] << std::endl;
        os << "}" << std::endl;
        return os;
    }
};

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * @brief 线性求解器
 *
 */
class LinearSolver
{
  public:
    /**
     * @brief 是否按照单位上三角阵求解.
     * @return
     * @param A
     * @param b
     * @param unit
     * @see
     */
    Vector SolveUpper(const Matrix &A, const Vector &b, bool unit = false);

    /**
     * @brief 是否按照单位下三角阵求解.
     * @return
     * @param A
     * @param b
     * @param unit
     * @see
     */
    Vector SolveLower(const Matrix &A, const Vector &b, bool unit = false);

    /**
     * @brief 是否选主元求解.
     * @return
     * @param A
     * @param b
     * @param choose
     * @see
     */
    Vector Solve(const Matrix &A, const Vector &b, bool choose = false);

    /**
     * @brief 是否优化为 LDL 分解.
     * @return
     * @param A
     * @param b
     * @param improve
     * @see
     */
    Vector Cholesky(const Matrix &A, const Vector &b, bool improve = false);

    /**
     * @brief Jacobi 迭代法.
     * @return
     * @param A
     * @param b
     * @param omega
     * @see
     */
    Vector Jacobi(const Matrix &A, const Vector &b, double omega = 1);

    /**
     * @brief Gauss-Sediel 迭代.
     * @return
     * @param A
     * @param b
     * @see
     */
    Vector GaussSediel(const Matrix &A, const Vector &b);

    /**
     * @brief 超松弛迭代.
     * @return
     * @param A
     * @param b
     * @param w
     * @see
     */
    Vector SOR(const Matrix &A, const Vector &b, double w);

    /**
     * @brief QR 分解法.
     * @return
     * @param A
     * @param b
     * @see
     */
    Vector QR(const Matrix &A, const Vector &b);

    /**
     * @brief 共轭梯度法.
     * @param A
     * @param b
     * @param x
     * @see
     */
    void CG(const Matrix &A, const Vector &b, Vector &x);

    /**
     * @brief 幂法。计算 A-muI 的最大特征值，然后加回 mu.
     * @return
     * @param A
     * @param mu
     * @see
     */
    double Power(const Matrix &A, double mu = 0);

    /**
     * @brief 反幂法。计算指定特征值的特征向量（半次迭代法）.
     * @return
     * @param A
     * @param lambda
     * @see
     */
    Vector InversePower(const Matrix &A, double lambda);

    /**
     * @brief 隐式 QR 方法.
     * @return
     * @param A
     * @param Q
     * @param ifSym
     * @see
     */
    Matrix ImplicitQR(const Matrix &A, Matrix *Q = nullptr, bool ifSym = false);

    /**
     * @brief 过关 Jacobi 方法。将 A 变换为特征值对角阵，返回正交阵
     * Q，每一列是特征向量，sigma > n 为关值减小的比例.
     * @return
     * @param A
     * @param sigma
     * @see
     */
    Matrix PassJacobi(Matrix &A, double sigma);

    /**
     * @brief 二分法。计算三对角对称阵的全部特征值.
     * @return
     * @param A
     * @see
     */
    Vector Bisection(const Matrix &A);

    /**
     * @brief 奇异值分解。返回分解后的对角阵，正交阵存放在 U,V 中.
     * @return
     * @param A
     * @param U
     * @param V
     * @see
     */
    Matrix SVD(const Matrix &A, Matrix *U = nullptr, Matrix *V = nullptr);
};

} // namespace CAGD
